Glassy dynamics in granular compaction: sand on random graphs 
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We discuss the use of a ferromagnetic spin model on a random graph to model granular com- 
paction. A multi-spin interaction is used to capture the competition between local and global 
satisfaction of constraints characteristic for geometric frustration. We define an athermal dynamics 
designed to model repeated taps of a given strength. Amplitude cycling and the effect of perma- 
nently constraining a subset of the spins at a given amplitude is discussed. Finally we check the 
validity of Edwards' hypothesis for the athermal tapping dynamics. 



I. INTRODUCTION 

Granular matter and glasses share a number of properties - such as off-equilibrium dynamics, aging, and hysteresis 
- and analogies between them have long jjj been pointed out. However it was not until the seminal experiments of 
the Chicago group on granular compaction were carried out that serious attempts were made to quantify such 
analogies. The experiments focus on the compaction behaviour of a large number of grains subject to repeated tapping 
and have become a paradigm for subsequent theoretical models. 

These models fall into roughly two classes: lattice-based models (|-||] in a finite-dimensional space (which in general 
do not admit analytic solutions) or mean field models ||j7| (where each site interacts with a large number of other 
sites). In this paper we discuss how models on random graphs may be used to describe aspects of the behaviour of 
granular matter which depend on the finite connectivity of the (disordered) grains, while still remaining analytically 
accessible. 

The aim of this paper is twofold: We discuss and motivate a simple spin-model defined on a random graph introduced 
as a model of granular compaction in ^| . In this model the random close packing density reached asymptotically after 
a large number of taps is identified with a dynamic phase transition. Secondly, we discuss an athermal dynamics ^,|), 
consisting of alternating periods of thermal dynamics at a certain temperature and quenches at zero temperature 
which take the system to a metastable state. 

Next, the compaction curve of the tapping process will be discussed and its two main features, the single particle 
relaxation threshold and the random close packing density (dynamical transition), will be analyzed in detail. We 
then present results of numerical simulations of tapping-amplitude cycling with a discussion of hysteresis and the 
asymptotic state of the model. Finally we investigate the statistical mechanics of the blocked configurations and 
discuss the validity of the so-called Edwards measure || . 

II. RANDOM GRAPH MODELS AND GRANULAR MEDIA 

A random graph jl0| consists of a set of nodes and bonds, with the bonds connecting each node at random to a 
finite number of others, thus from the point of view of connectivity appearing like a finite-dimensional structure. Each 
bond may link up two sites (a graph) or more (a so-called hypergraph). 

Formally, a random graph of N nodes and average connectivity c is constructed by considering all N(N — l)/2 
possible bonds between the nodes and placing a bond on each of them with probability c/N. In other words the 
connectivity matrix is sparse and has entries 1 (bond present) and (no bond), which are independent and 
identically distributed variables with probability c/N and 1 — c/N respectively. The resulting distribution of local 
connectivities is Poissonian with mean and variance c. The resulting structure is locally tree-like but has loops of 
length of order ln(A). Although there is no geometric concept of distance (in a finite dimensional space), a chemical 
distance may be defined by determining the minimum number of steps it takes to go from one given point to another. 

In a similar fashion, graphs - strictly speaking hypergraphs - with plaquettes connecting 3 or more nodes each 
may be constructed. Choosing Cy^ = 1(0) randomly with probability 2c/ N 2 (1 — 2c/ N 2 ) results in a random 3- 
hypergraph, where the number of plaquettes connected to a site is distributed with a Poisson distribution of average 
c. An illustration of part of such a graph is shown in Figure ^. 



1 



FIG. 1. A part of a random graph (strictly speaking a hypergraph) with triplets of sites forming plaquettes illustrating its 
locally tree-like nature (no planarity or geometric sense of distance are implied). 



Spin models on random graphs have been investigated for almost 20 years [ pd| since they may be considered as being 
halfway between infinite-connectivity models and finite-dimensional models, having to a certain extent the analytic 
accessibility of the former within the framework of mean-field theory, yet the finite connectivity of the latter. Interest 
in these models has intensified lately since they occur in the context of random combinatorial optimization problems 
jl2| and inroads have been made towards their analytic treatment beyond replica-symmetry |L%|l4| ]. 

In the context of modelling the compaction of granular matter, random graphs are the simplest structures with a 
finite number of neighbours. The finite connectivity is a key property, which goes beyond the simple fact that the 
grains in a granulate are in contact with a finite number of neighbouring grains. For instance kinetic constraints, which 
are a prominent feature in many models of granular behaviour |l5| can only be meaningfully defined in models with a 
finite connectivity. Furthermore cascades found experimentally during the compaction process may be explained by 
interactions between a finite number of neighbouring sites, where one local rearrangement sets off another one in its 
neighbourhood and so on ||. 

Another reason for the use of random graphs lies in the disordered structure of granular matter even at high 
densities. A random graph is the simplest object where a neighbourhood of each site may be defined, but has no 
global symmetries like a regular lattice. Additionally, the locally fluctuating connectivity may be thought of as 
modelling the range of coordination numbers of the grains DM . 



III. THE MODEL 



In the following we consider a 3-spin Hamiltonian on a random hypergraph where N binary spins Si = ±1 interact 
in triplets 

H=-pN = - CtjkStSjSk (1) 

i<j<k 

where the variable Cy^ = 1 with i < j < k denotes the presence of a plaquette connecting sites i, j, k and = 
denotes its absence. 

This Hamiltonian has recently been studied on a random graph in the context of satisfiability problems in combi- 
natorial optimization |l7| , on a random graph of fixed connectivity |l4| , and on a 2D triangular lattice [ 
a trivial ground state where all spins point up and all plaquettes are in the configuration + + + giving a contribution 

of —1 to the energy. Yet, locally, plaquettes of the type h, — I — , H (satisfied plaquettes) also give the same 

contribution. This results in a competition between local and global satisfaction of the plaquettes: Locally any of the 
satisfied plaquettes are equivalent (thus favouring a paramagnetic state), yet globally a ferromagnetic state may be 

favoured, since there are few configurations satisfying all plaquettes, where 4 configurations + + + , h, — + — , H 

occur in equal proportions. For c > c c ~ 2.75 |17[, ground states have a positive magnetisation, which may be inter- 
preted as the onset long-range order and of a possibly crystalline pcjplj state of the granular medium. 
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FIG. 2. The phase space of three spins connected by a single plaquette. Configurations of energy —1 (the plaquette is 
satisfied) are indicated by a black dot, those of energy +1 (the plaquette is unsatisfied) are indicated by a white dot. 

However two spin flips are required to take a given plaquette from one satisfied configuration to another. Thus 
an energy barrier has to be crossed in any intermediate step between two satisfied configurations as illustrated in 
figure ^. In the context of granular matter this mechanism aims to model the situation where compaction follows a 
temporary dilation; for example, a grain could form an unstable ('loose') bridge with other grains before it collapses 
into an available void beneath the latter. This mechanism, by which an energy barrier has to be crossed in going 
from one metastable state to another, has recently been argued to be an important ingredient in models of granular 
compaction pl[ . 

The crucial feature of the model responsible for the slow dynamics, however, is the degeneracy of the four config- 
urations of plaquettes with SiSjSk — 1 resulting in competition between satisfying plaquettes locally and globally. In 
the former case, all states with even parity may be used, resulting in a large entropy and in the latter, only the + + + 
state may be used. A dynamics based on local quantities will thus fail to find the magnetised configurations of low 
energy. 

This mechanism has a suggestive analogy in the concept of geometrical frustration of granular matter, if we think 
of plaquettes as granular clusters. When grains are shaken, they rearrange locally, but locally dense configurations 
can be mutually incompatible. Voids may appear between densely packed clusters due to mutually incompatible grain 
orientations between neighbouring clusters. The process of compaction in granular media consists of a competition 
between the compaction of local clusters and the minimisation of voids globally. 



A. Modelling tapping 

There have been many kinds of dynamical schemes to model the behaviour of granular media under tapping. A 
recurring theme is the alternation of periods of randomly perturbing the system and periods in which the system 
is allowed to settle into a mechanically stable state. These have included nonsequential Monte Carlo reorganisation 
schemes , the ratio of upward to downward mobility of particles on a lattice Q , or variable rates of absorption 
and desorption p3| . 

In the same spirit, we treat each tap as consisting of two phases. First, during the dilation phase, particles are 
accelerated and are relatively free to move with respect to each other for a time. In the second phase, the quench 
phase, particles relax until a mechanically stable configuration is reached. 

As initial condition we use a configuration obtained by quenching the system from a configuration where the spins 
are chosen independently to be ±1 with equal probabilities. To mimic the action of tapping, we choose the following 
dynamics of the spins. The dilation phase is modelled by a single sequential Monte-Carlo-sweep of the system at a 
dimcnsionless temperature T: A site i is chosen at random and flipped with probability 1 if its spin si is antiparallcl 
to its local field hi, with probability exp(— hi/T) if it is not, and with probability 0.5 if hi = 0. This procedure is 
repeated N times. Sites with a large absolute value of the local field h{ thus have a low probability of flipping into the 
direction against the field. Such spins may be thought of as being highly constrained by their neighbours. This differs 
somewhat from the dilation-phase model in I^Q, where a certain fraction of spins is flipped regardless of the value 
of their local field. We claim that our present dynamics is rather more realistic in the context of vibrated granular 
media; if grains are densely packed (strongly 'bonded' to their neighbours), they are less likely to be displaced during 
the dilation phase of vibration than grains which are loosely packed. 

The quench phase is modelled by a quench of the system at T — 0, which lasts until the system has reached a 
blocked configuration, i.e. each site i has Si = sgn(/ii) or hi = 0. Thus at the end of each tap the system will be in a 
blocked configuration. 
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This dynamics may be thought of as a series of quenches where the initial condition for each quench is obtained by 
perturbing the result of the previous quench. It is a simplified version, suitable for spin models, of the tapping dynamics 
used in cooperative Monte Carlo simulations of sphere shaking p2fl . In the context of combinatorial optimization it 
corresponds to the class of random-restart algorithms (e.g. PjJ ), which include some of the most efficient algorithms 
for the solution of optimization problems. 

IV. THE COMPACTION CURVE 

An example of a single run of the system is shown in figure |3|. 

We can identify three regimes of the dynamics: first, a very fast increase of the density up to a density po during 
the first tap, then a slow compaction regime which takes the density up to poo, and finally an asymptotic regime. 

In the first regime, all sites orient their spins in parallel with the local field acting on that site. This quench 
corresponds to a fast dynamics whereby single particles locally find the orientation maximizing the density leading to 
the density of po [p5[ . In Q this density was termed the single-particle relaxation threshold (SPRT). 




10° 10 2 10 4 10 6 

taps 

FIG. 3. Compaction curve at connectivity c = 3 for a system of 10 4 spins with T = .4. The data stem from a single run and 
the fit (smooth solid line line) follows (^) with parameters poc = .989, po = .843, D = 4.716, and r = 52.46. The long-dashed 
line (top) indicates the approximate density 0.954 at which the dynamical transition occurs, the long-dashed line (bottom) 
indicates the approximate density 0.835 at which the fast dynamics stops, the single-particle relaxation threshold. 

The second phase of the dynamics consists of removing some of the remaining frustrated plaquettes and gives a 
logarithmically slow compaction Pj5l leading from a density po to poo. The resulting compaction curve may be fitted 
to the well-known logarithmic law H] 

p(t) = Poo - (Poo - Po)/(l + l/D ln(l + t/r)) , (2) 

which may also be written in the simple form 1 + t(p)/r = cxp {D £ } , implying that the dynamics becomes 

slow (logarithmic) as soon as the density reaches po- In this regime, most spins have a nonzero local field acting on 
them, which keeps them fixed in a certain direction [ p6[ . The corresponding grains are firmly held in place by their 
neighbours. However, during the dilation phase some of them have their orientation altered, altering the local fields 
acting on their neighbouring grains by a finite amount, which could cause them to flip in turn. The dynamics of the 
grains with zero local field may alter the local field of their neighbours, and induce a previously blocked grain to flip. 
In this way a cascade || of flips may ensue. 

With increasing density, free-energy barriers rise up causing the dynamics to slow down according to (^|) . The point 
where the height of these barriers scales with the system size marks a breaking of the ergodicity of the dynamics, a 
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break-up of the phase-space into a large number (scaling exponentially with the system size) of disconnected clusters, 
and a saturation of the compaction curve. For small driving amplitudes, we thus identify the asymptotic density 
(random close packing) with a dynamical phase transition |t],|27]-^] . 

In the following we will examine in detail the SPRT and the dynamical transition. 

A. The single-particle relaxation threshold 

The first tap is modelled as a quench at zero temperature. At the end of this, each site is connected to more (or as 
many) unfrustrated plaquettes than frustrated ones. The spin of any site where this is not the case would flip under 
the zero-temperature dynamics, turning frustrated plaquettes into unfrustrated ones. The question of the density 
reached after a quench from random starting conditions is highly non-trivial, since its resolution involves the basins 
of attraction of the zero-temperature dynamics. 

The problem may be illustrated by considering a single site i connected to 2ki other sites and subject to the local 
field hi = l/2^jjt CijkSjSk- For random initial conditions, the values of h — hiSi are binomially distributed with a 
probability of Cfe _ ; s, 2 (l/2) ki if fcj — ij is even and zero if it is odd. If U < zero-temperature dynamics will flip this 
spin, turn li to —U and turn (hi ± h)/2 satisfied (dissatisfied) plaquettes connected to it into dissatisfied (satisfied) 
ones. This will cause the lj of ki ± li neighbouring sites to decrease (increase) by 2. This dynamics stops when all 
sites have I > 0, giving p — 1/(3 A) J2i h- 

This process is made complicated by correlations between the local fields of neighbouring sites. Neglecting these 
correlations we arrive at a simple population model of A units, each with a Poisson distributed value of ki and a 
value of h distributed according to the initial binomial distribution. At each step a randomly chosen element with 
negative h has its h inverted, and ki ± U randomly chosen elements have their values of I decreased (increased) by 
2 until U > Vi. This simplistic model works surprisingly well at low values of the connectivity c (with an error of 
about 10% up to c = 6), but obviously fails completely at large values of c or in fully connected models. 

In principle the differential equations describing the population dynamics could be solved analytically. Here we 
simply report the results for running the population dynamics numerically with A = 10 4 at c = 3. We obtain 
Po = 0.835(1) which is shown as a dotted line in Figure [| Note that this density is found to be much higher than that 



of a typical 'blocked' configuration with k > Vz which is found to be 0.49 (see section VI and also the discussion in 
p(i||). Despite the fact that these 'blocked' configurations are exponentially dominant, the total basin of attraction of 
the configurations at po dominates the space of random initial conditions. 

Another significant feature of this regime is that a fraction of spins is left with local fields exactly equal to zero, 
which thus keep changing orientation p6| . These spins may be compared to so-called rattlers, pi]] , i.e. grains which 
change their orientation within well-defined clusters []l6) . These will be used as a tool to probe the statistics of blocked 



configurations in section VI 



To conclude this section, the SPRT density appears as the density which is reached dynamically by putting each 
particle into its locally optimal configuration, as has also been found in lattice-based models and simulations of 
sphere packings (|l],^Cj], which show both fast and slow dynamics. 

B. The dynamic transition 

The dynamical transition is marked by the appearance of an exponential number of valleys in the free-energy land- 
scape and thus a breaking of ergodicity p7|-p9|] . In the event that the dynamics is thermal, equilibration times diverge 
at the temperature corresponding to the dynamic transition. Cooling the system down gradually from high tempera- 
tures will also result in the system falling out of equilibrium at the dynamical transition temperature. Furthermore, 
the energy will get stuck at the energy at which the transition occurs. 

Since this phenomenon is the result of the drastic chang e in th e geometry of phase space, it is not surprising that 



we also find it in the athermal dynamics defined in section [II A . Either gradually decreasing the tapping amplitude 
T or tapping at a low amplitude for a long time will get the system to approach the density (energy) at which the 
dynamical transition occurs. We thus identify the random close packing limit in this model with a dynamical phase 
transition. 

To support this picture, we give a simple approximation for the density poo at which the dynamical transition 
occurs. Using the replica-trick \nZ = lim„^o d n Z n |3^] and standard manipulations, we obtain for the average of the 
n-th power of the partition-function of the Hamiltonian (mj averaged over the ensemble of random graphs 
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((Z n ))=J[j\(a)exp^-N fc>(Wc(*)) ( 3 ) 
+c/3-c/3 ^ c{uj)c{f)c(a)exp{P^2uj a T a a a } 

where c(a) is an order parameter function defined on the domain of the 2 n vectors a a — ±1. 

The general replica-symmetric ansatz is incorporated in c(a) = J dhP(h) ( 2 cosii(/3H))" • Taking P(/i) = <5(/i) gives 
the paramagnetic solution, valid in the high-temperature phase, resulting in a free energy 

/3/(/3) = -c/31n(cosh/3) + ln(2) . (4) 

To determine the temperature at which the dynamics transition occurs, a replica-symmetry breaking (RSB) ansatz 
is required (p^,p8l. A simple variational ansatz fl33"|,p"7 34 implementing one step of RSB given by 



„ _ ^ f J dh b G A (h b )e" hb ^-^ y 

a) 11 l / dh b G A (h b )[2cosh(f3h b )] m f ' ' 



n/rti 

c\n ) = 

6=1 



where G A (h) is a Gaussian with zero mean and variance A, gives the free energy subject to the variational ansatz 
f((3) = extr A ,m/i(/3, A,m) with 

aclQK , |^z(/3\/Az)[2cosh(/3VAz)] m - 1 sinh(/3yAz) 

ph(p,A,m) = ; p= 

V ' /L»z[2cosh(/VAz)]™ 

-— - ln( J L>z[2cosh(/3\/Az)] m ) - c/(3m) ln( / J J Dz l Dz 2 Dz i 

8cosh(/3v / Azi) cosh(/3v / Az 2 ) cosh(/3v / Az 3 ) cosh(/3) + 8 sinh(/3\/Azi) sinh(/3v / Az 2 ) sinh(/3\/Az3) sinh(/3) 



where -D(^) denotes the Gaussian measure with zero mean and variance one. The dynamical transition occurs at 
a temperature |3^] where d(/3f(j3,A,m))/dm evaluated at m = 1 develops a minimum at finite A [ p7p^| . The 
corresponding density is marked with a horizontal line in Figure |3| and agrees well with the asymptotic density 
reached by the tapping dynamics. However this asympotitic density is not the highest density, that can be reached 
without putting the system into an ordered configuration, however it is the highest such density which is reached by 
a local dynamics. 



V. AMPLITUDE CYCLING AND THE STATIONARY STATE 



The tapping dynamics introduced in section III A may be used to increase and decrease the tapping amplitude 
successively. This amplitude cycling is an important protocol in real and numerical experiments. The ramp rate 
is defined as the ratio 5 (Amplitude) / t where r is the number of taps spent at each amplitude, which is changed with 
an increment of 6 (Amplitude) after each series of taps. 

The results of increasing T continuously from to 2 and back again at two rates 10 -4 and 10 -5 per tap is shown 
in figure |J; here, clearly t = 1, and the ramp rate is in fact just S(Amplitude). As expected, both at high and at 
low cycling rates, the density first reaches the SPRT po, then increases with increasing amplitude and time, until it 
decreases again at large values of T. As the amplitude is decreased, the system reaches p^. The part of the curve 
where T is increased for the first time is conventionally || called the irreversible branch, while the reversible branch 
refers to the section where T is subsequently decreased and then increased again. 



However the similarity of this simple picture with experimental results of |2| is deceptive. From section IV B we 
know that at fixed, finite but low amplitudes the model eventually reaches poo. As a result, the branches of increasing 
amplitude at low T do not coincide for high and low rates of change of the amplitude. At low rates of change the 
density as a function of T is higher than at high rates of change. This irreversibility of the so-called reversible branch 
is thus due to the system not reaching a steady-state at each value of the amplitude. This behaviour has also been 
observed in other models P6LH . 
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FIG. 4. The results of ramping the amplitude up and down again at two different rates, 10 -4 (solid line) per tap and 10~ 5 
(dashed line). The lower rate results in a steeper increase of the density during the increase of the amplitude (lower branches 
of solid and dashed lines). 



In fact, in the limit of infinitely slow increments of T, the irreversible branch would disappear, and p would become 
a single- valued function of T. This would be in direct contradiction to the experimental results of 0], where at 
low amplitudes the steady-state density cannot be reached with a sufficiently large number of taps, at least within 
experimentally realisable times. 

This phenomenon may be thought of as follows: Some particles in a granular assembly are so strongly constrained 
that they will never (at least within experimentally realisable times) be moved by taps of a sufficiently low amplitude. 
In the dynamics of our model described thus far, however, sites with a high local field may be flipped at any finite 
value of T with a correspondingly small but finite probability, leading the system eventually to . 

To model this effect, it is not sufficient to prevent spins with a large value of the local magnetic field from being 
flipped since the dynamics of their neighbours will eventually lead to a reduction of their field, freeing the previously 
constrained spins. Instead, we assign to each site i a real number n between zero and one, and during the dilation 
phase of the dynamics, we only flip spins at sites i with rj < T. During the 'quench' phase any spin may be flipped. 

In the language of grains, rj represents the strength with which a grain is constrained by its neighbours; sites i such 
that Ti > T will be permanently resistant to being displaced at an intensity of vibration T. In a real system, these 
thresholds would be determined by details of the inter-grain force netw ork. 

The aim of this modification is to check if the scenario of section IV B survives, since in principle it is possible that a 
new and lower value of p^ emerges after the amplitude has been increased and decreased: At high shaking amplitudes 
frustrated plaquettes (clusters) might be generated, which cannot be eliminated at low values of the amplitude, similar 
to the scenario proposed in |37| . 



7 




Figure g, where the cycling of Figure |4| is repeated with the constraining of spins, shows that p^ remains unaltered. 
In fact, Figures ^ and look remarkably similar. The effect of constraining some spins at low values of T emerges 
when the ramp rate is decreased substantially, in particular by allowing the spins to 'equilibrate' at each amplitude 
by choosing a large r. In figure^ we thus increase the shaking amplitude T from 0.2 to 2 in steps of 0.2 with t = 10 7 
taps waited at each amplitude step. This makes sure that a steady-state of the density has been reached at each 
amplitude. We find, in this case, that at low amplitudes the immobile spins cause the system to reach a steady-state 
with a density lower than p^. Despite the fact that in p6|,p[, very low ramp rates were used, with large 'waiting 
times' r at each tap, this behaviour was not observed; rather the results of all these simulations implied that the 
asymptotic density p,^ would always be approached in the limit of sufficiently low ramp rates. 

One may view the (random) configuration of the immobile spins at each value of T as an additional quenched 
disorder and their effect on neighbouring mobile spins, as a random local field. Presumably the dynamics in the sub- 
space of phase-space corresponding to the mobile spins (with fixed local fields due to the immobile spins) undergoes a 
dynamic transition as the corresponding steady-state density is reached. The result that p^ is reached decreasing T 
from above even though a finite fraction of spins has been rendered immobile at low temperatures, is quite remarkable: 
it is a testament to the paramagnetic nature of the model at densities below poo . In a glassy state, one would expect 
the configurations of spins reached at high values of T and subsequently frozen to alter the behaviour of the system at 
lower values of T. The paramagnetic system manages very well to adapt the mobile spins to the configuration of the 
immobile ones produced at higher values of T - not however to random configurations of the immobile spins, which 
are responsible for the irreversible branch. 

These results demonstrate a rather fundamental difference between thermal excitations in glassy systems and inten- 
sities of mechanical vibration in granular media. In the glassy phase of a system, one would expect the configurations 
of spins reached at high values of temperature and subsequently frozen, to alter the behaviour of the system at lower 
values of temperature. In granular media, however, it is important to let the system reach the asymptotic density at 
each value of the shaking intensity! 1 , in order even to begin to observe the hysteresis that must result when mobile 
grains become part of immobile clusters |3q| , generating a type of 'quenched' disorder at least at low vibrational 
intensities. The difference between Figs. and ^ clearly illustrates this. Given that the experiments results 
were done in such a way that the system was allowed to reach the asymptotic density at each value of the tapping 
amplitude, our results indicate that 'jamming' |}9) of grains caused by the force network might be responsible for the 
fact that the p^ is not reached by tapping solely at low amplitudes. 
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VI. BLOCKED CONFIGURATIONS AND THE EDWARDS-MEASURE 



In this section we focus on the statistical mechanics of the blocked configurations referred to earlier, and use these 



results to address the question of ergodicity of the tapping dynamics. After each tap according to section [II A , 
the system is in a blocked configuration, i.e. each site i has s$ = sign(/ii) or hi = 0. The Edwards hypothesis ^C] 
states that in the steady state along the reversible branch all mechanically s table configurations at a given density 



are equiprobable. We test this hypothesis for the tapping dynamics of section [III A 

We begin by calculating the average entropy of blocked configurations at a given density. In principle, we would 
need to average the logarithm of the number of blocked states over the ensemble of random graphs; this so-called 
quenched average can be expected to be self-averaging. For simplicity, we restrict ourselves to the so-called annealed 
average and compute 

s annealed(p) = Jj ln ((-^(p)>> > ^quenched (p) = ^(( ln W(fi)))) > (6) 

which gives an upper bound to the quenched average. The number of blocked configurations M{p) may be written 
easily as 

AA (/ 5 ) = Il( £ £ 5(h i ;l/2j2Ci j kS j s k )e(h i s i ))s[p-l/(3N)J2 h i s i) > ( 7 ) 

, Si— ±1 hi — — oo j.k 

where S(x; y) = 1 if x = y and otherwise, denotes a Kronecker-delta and Q(x) denotes a discrete Heaviside step 
function with Q(x) = 1 if x > and otherwise. After using integral representations for the Kronecker-deltas and 
standard manipulations, one easily obtains 



blocked / \ _ „„ f „ 
s armealedW - extr a , &i/3 



-pp + 8c/3 (a 3 + b 3 ) - c/3 + In 2 ^{e^ 3 a/b) h I h {Acab) + 2I (Acab) 



h=l 



(8) 



where Ih(x) denotes the modified Bessel- function of the first kind of order h. In Figure |^, the entropy of blocked 
configurations s a nnealed(/ 9 ) ^ s sri0wn along with the paramagnetic entropy given by ( Jll| ) below, which is derived by 
including terms with negative Sihi- 
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FIG. 7. The paramagnetic entropy (top) and the entropy of blocked states in the annealed approximation (bottom) for c = 3 
as a function of the density p. The data with errorbars shows the results of exhaustive enumerations of a system with N — 28 
averaged over 100 samples. The results for the paramagnetic state are marked with squares. 

These expressions were evaluated for c — 3 and the results are shown in figure 0. From this data, the lowest density 
at which blocked configurations occur (s annea i e{ j(/o) = 0) is found to be p = .09, and the density of a randomly chosen 
blocked configuration (the maximum of s a nnealed(/°)) is P = -49. 

Similarly, one may calculate the fraction g of connected sites with zero local magnetic field in a blocked configuration. 
As we will argue below, this is a useful quantity to test the Edwards hypothesis. ^From (|h one obtains 
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I (4cab) 



EZi(e fj/3 a/b) h h(4cab) + I (4cab) 



(9) 



where for any given value of p, the values of (3 and of the order parameters a and b are given by the extremisation 
condition in equation (|J), and where we have subtracted a trivial term e~ c corresponding to the fraction of unconnected 
sites. 

Analogously, the fraction of connected sites with zero local magnetic field at a given density (without the blocking 
condition) is given by 



9' = 



h(Acab) 



Eh=i [{e 0/3 a/b) h + {e-P/%/a) h ] I h {Acab) + I {4,cab) 
where the values of a, 6, and/3 follow from extremizing over 



(10) 



s annealed (^) — ex ^ T a,b,0 



-I3p + 8c/3 (a 3 



/ oo -\ 

c/3 + In j 2 [(e 0/3 a/b) h + (e^fe/a)' 1 ] 4(4ca6) + 2I (4cab) 



(11) 



Figure |^ shows g and gi as a function of p. As expected, both quantities decrease monotonically with p. Again, we 
also give the results of exhaustive enumerations of a system with N = 28, averaged over 100 samples in order to test 
the validity of the annealed approximation. The annealed curve for blocked states and the numerical results show 
significant differences. For this reason, we also give the result of the much more involved replica-calculation of the 
quenched average JO] as the dashed line, which agrees very well with the numerical results. 
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FIG. 8. The fraction of connected spins with zero magnetic field plotted against p for c = 3. The top solid line gives the 
analytical result g (||) for blocked configurations, the bottom solid line gives that for the paramagnetic state gl (jn]) (without 
the blocking condition). The results of the quenched average are shown as a dashed line. The data with errorbars shows the 
results of exhaustive enumerations of a system with N = 28 averaged over 100 samples. The results for the paramagnetic state 
are marked with a square. 



With these results, we may now test the Edwards hypothesis for this model under the tapping dynamics of section 



[II A . If blocked states at the asymptotic density are accessed with equal probability a plot of the fraction of connected 
sites with zero local field versus the density should coincide with the results of (|^) at the asymp totic density. 

The value of g is a useful quantity to test this hypothesis, since as discussed in section IV A, the fraction of spins 
with a given local field is intimately related to the fast quenching dynamics. Also, as a one-time quantity it may be 
measured easily. 

In Figure ||, we show the results of 4 single runs at T = 0.4,0.56,0.7,1.5, plotting g against p. It shows clearly 
that at the asymptotic density, the (quenched) result for g against p and the results of the tapping dynamics agree to 
within numerical accuracy of the analytical result, Wc tentatively conclude that the Edwards hypothesis is valid in 
this model at low tapping intensities and at the asymptotic density reached by low-intensity tapping. Further results 
will be reported elsewhere [[ll| . 

During the compaction phase, however, the blocked configurations accessed dynamically have a lower value than that 
of the exponentially dominant blocked configurations contributing to (^|). The result that the blocked configurations 
accessed dynamically have a lower value of g than the typical ones at the same value of p, may be explained to a certain 
extent as follows. Configurations with a large value of g are favoured cntropically, since spins with zero magnetic field 
may be flipped leaving the configuration blocked (provided this does not cause the local field of their neighbours to 
change sign). In the quenching dynamics, however, there is no such mechanism; sites with zero magnetic field are 
created only by spin-flips of neighbouring sites. 

Nevertheless it is remarkable that for the three lower values of T (light tapping), where the asymptotic density is 
very close to poo, the three traces nearly fall onto a single line, indicating that also during compaction the blocked 
configurations are sampled according to a certain ensemble which depends on the density only. This occurs even 
though the plots of rho versus the number of taps do not coincide for these amplitudes. 
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FIG. 9. The fraction g of connected sites with zero local magnetic field during 4 runs of a system with N = 10000, c = 
at T = 0.4 (dots), T = 0.56 (+), T = 0.7 (x), and T = 1.5 (circles). The solid line gives the analytic annealed result of (§) 
the dashed line the corresponding quenched result. The lines left and right indicate the approximate values for po and p a 
respectively. 



VII. CONCLUSION 



To conclude, we have presented a finitely connected spin model of vibrated granular matter, where we have built 
upon earlier work We argue that spin-models on random graphs may serve as models of granular matter, since 
they show no symmetries in the way a regular lattice does. They also arise as the Bcthe-lattice approximation to finite- 
dimensional models. Multi-spin interactions generically arise when models of geometric frustrations are transferred 
to the Bethe-lattice. We discuss one of the simplest models of this class, the ferromagnetic 3-spin model. Due to 
competition between satisfying the interactions globally and locally the model never reaches the ferromagnetic state. 
This mechanism aims to model the geometric frustration incurred by packings which arise from maximising the density 
locally. 

We also discuss the problems incurred by glassy models in the context of amplitude cycling. In order to model the 
effect of some grains being rendered completely immobile by large intergranular forces, we investigate the effect of 
constrainin g 'b y hand' some of the spins in the context of amplitude cycling experiments. We also test the Edwards 
hypothesis [|f0|| in the context of our model, and concluded that it is valid for the tapping dynamics used. 

Acknowledgements: It is a pleasure to thank S. Franz, F. Ricci-Tersenghi, and R. Zecchinafor 
illuminating discussions. 
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